home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / msecant.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-27  |  2.7 KB  |  109 lines
  1. /*
  2. Compute periodic orbit of a period r by secant method
  3. (Solve F(x)=f^r(x)-x-xoffset= 0)
  4. ------------------------------------------------------
  5. Note:    x[0...n-1] convention 
  6.     All computation in Euclidean coords
  7.     
  8. Input:    x[n]= current guess, n= phase space dimension, r= period
  9.     tolx= tolerated error in x (sum), tolf= tolerated error in F(x) (sum)
  10.     ntrial = # of maximum secant step
  11. Output:    x[n]=new guessed solution, *ndid= acutal # of secant step done,
  12.     *err= actual error in x (sum)
  13. */
  14.  
  15. #define DELTAX 1.e-6
  16. #define FREERETURN {*ndid = k;free_dmatrix(alpha,0,n-1,0,n-1);free_dvector(x2,0,n-1);free_dvector(xt,0,n-1);free_dvector(beta,0,n-1);\
  17.     free_ivector(indx,0,n-1);return;}
  18.  
  19. void msecant(x,ntrial,n,r,tolx,tolf,err,ndid)
  20. int ntrial,n,r,*ndid;
  21. double x[],tolx,tolf,*err;
  22. {
  23.     int i,j,k,*indx,*ivector();
  24.     double tolx10,fabs(),errf,d,*x2,*xt,*beta,**alpha,*dvector(),**dmatrix();
  25.     void usrfun(),usrfun2(),usrfun3(),ludcmp(),lubksb(),free_dmatrix(),free_ivector(),free_dvector();
  26.     extern int stop,enable_period,cur_color,model,mapping_on,fderiv_on;
  27.  
  28.     indx = ivector(0,n-1);
  29.     x2 = dvector(0,n-1);
  30.     xt = dvector(0,n-1);
  31.     beta = dvector(0,n-1);
  32.     alpha = dmatrix(0,n-1,0,n-1);
  33.  
  34.     tolx10 = tolx /10;
  35.     *err = 1.e30;
  36.     /* initial second guess */
  37.     for(i=0;i<n;i++) x2[i] = x[i] + DELTAX;
  38.  
  39.     for(k=0;k<ntrial;k++){
  40.         if(mapping_on){
  41.             if(fderiv_on)
  42.                 usrfun(x,alpha,beta,r,n);
  43.             else
  44.                 usrfun2(x,x2,alpha,beta,r,n);
  45.         }
  46.         else {
  47.             usrfun3(x,x2,alpha,beta,r,n);
  48.         }
  49.         if(enable_period==1){
  50.             dist_periodic(xt,beta,n);
  51.             for(i=0;i<n;i++)beta[i]=xt[i];
  52.         }
  53.         /*
  54.         printf("NTRIAL=%d\n",k);
  55.         printf("usrfun x: ");
  56.         for(i=0;i<n;i++) printf("%g ",x[i]);    
  57.         printf("\n");
  58.         printf("usrfun f: ");
  59.         for(i=0;i<n;i++) printf("%g ",beta[i]);    
  60.         printf("\n");
  61.         printf("usrfun alpha:");
  62.         for(j=0;j<n;j++){
  63.             for(i=0;i<n;i++) printf("%g ",alpha[j][i]);    
  64.             printf("\n");
  65.         }
  66.         */
  67.  
  68.         errf = 0.0;
  69.         for(i=0;i<n;i++) errf += fabs(beta[i]);
  70.         if(errf <= tolf) FREERETURN
  71.  
  72.         ludcmp(alpha,n,indx,&d);
  73.         /* should be placed here not to proceed further
  74.         in case of singular matrix (singular matrix often
  75.         arise when all the machine accuracy was lost in
  76.         computing alpha. This happens in case of good
  77.         convergence,too) */
  78.         if(stop){
  79.             FREERETURN
  80.         }
  81.         /*
  82.         for(j=0;j<n;j++){
  83.             for(i=0;i<n;i++) printf("%g ",alpha[j][i]);    
  84.             printf("\n");
  85.         }
  86.         */
  87.         lubksb(alpha,n,indx,beta);
  88.  
  89.         /*
  90.         printf("new beta:");
  91.         for(i=0;i<n;i++) printf("%g ",beta[i]);    
  92.         printf("\n");
  93.         */
  94.         *err = 0.0;
  95.         for(i=0;i<n;i++){
  96.             *err += fabs(beta[i]);
  97.             x2[i] = x[i];
  98.             x[i] -= beta[i];
  99.         }
  100.         if(*err <= tolx) FREERETURN
  101.         for(i=0;i<n;i++) if(fabs(beta[i]) < tolx10) x[i] -= tolx10;
  102.  
  103.         /* draw intermediate steps */
  104.         (void) draw_record_pwf(x,cur_color,1,3,1,0);
  105.         (void) draw_record_pwf(x,cur_color,1,3,1,0);
  106.     }
  107.     FREERETURN
  108. }
  109.